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1.0  INTRODUCTION 


The  puipose  of  this  work  was  to  analyze  seismic-hydro¬ 
acoustic  signals  recorded,  at  an  ocean  bottom  geophone/ 
hydrophone  array,  during  a  shallow-water  explosive  shot 
experiment  conducted  by  the  Naval  Coastal  Systems  Center. 

This  analysis  was  to  employ  a  variety  of  special  signal 
processing  techniques,  in  order  to  determine  the  degree  of 
processing  gain  realizable  and  to  characterize  the  received 
signals  and  their  propagation  paths. 

In  the  experiment,  charges  having  the  equivalent  of  .8 
pounds  TNT  were  exploded  at  different  ranges  along  a  110  km 
southwesterly  radial  leg,  from  the  Naval  Coastal  Systems 
Center  Stage  I  array  located  in  the  Gulf  of  Mexico  near 
Panama  City,  Florida.  Signals  from  the  shots  were  received 
and  tape  recorded  at  the  array,  which  contained  seven  tri- 
axial  geophones  and  six  hydrophones  deployed  along  three 
radii.  Array  configuration  is  defined  in  Figure  la  and 
Table  1.  All  geophones  were  buried  about  .3  meters  deep  in 
the  ocean  bottom;  all  hydrophones  were  suspended  about  .6 
meters  above  the  bottom.  Water  depth  was  approximately  30 
meters  at  the  array,  gradually  increasing  to  about  40  meters 
at  range  25  kilometers.  Further  out,  the  bottom  slope  was 
about  .8  meters  downward  per  kilometer. 

Thirty  events  were  digitized,  representing  shot  ranges 
from  3.73  to  108.89  kilometers.  An  analysis  was  made  of 
these  digital  time  series,  using  reduced  travel  time-distance 
plots,  concentrating  on  the  high  signal -to-noise  events  within 
30  kilometers  of  the  array.  From  this  analysis,  a  preliminary 
seismic  propagation  model  was  established  to  explain  the 
observed  mode  structure.  A  subset  of  these  events  was  then 
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TABLE  1 

SENSOR  DISPLACEMENTS 


Sensor 

Distance  from  Array  Center  (meters) 

H1 

.635 

H2 

.457 

H3 

.762 

H4 

8.839 

H5 

8.966 

H6 

9.322 

G1  - 

6.045 

G2 

5.893 

G3 

6.121 

G4 

35.408 

G5 

34.925 

G6 

35.408 

G7 

0.000 

(From  Reference  1) 
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subjected  to  intensive  signal  processing,  including  time 
series  and  spectral  analysis,  coherent  crosscorrelation  and 
beamforming,  and  polarization  filtering. 


Substantial  signal - to-noise  ratio  (SNR)  gains  were 
achieved  by  use  of  the  coherent  crosscorrelation  techniques, 
although  certain  attributes  of  this  digitized  data  set  may 
be  unrepresentative  of  typically  encountered  signals.  In 
particular,  a  very  highly  correlated,  but  apparently  unphysi¬ 
cal,  signal  was  found  to  contaminate  the  data  set.  It  is 
not  clear  to  ENSCO  what  the  source  of  this  artifact  might 
have  been.  Application  of  polarization  filters  yielded 
generally  unusable  or  untenable  results,  even  though  the 
signals  clearly  should  have  been,  at  least  partially,  suscep¬ 
tible  to  this  processing.  A  possibility  exists  that  the 
transfer  of  code  and  modification  of  (previously  successful) 
poliarization  filter  algorithms,  for  the  present  application, 
was  not  made  satisfactorily. 

Both  seismometer  and  hydrophone  -  received  signals  were 
analyzed  as  part  of  this  investigation.  In  spite  of  the  fact 
that  the  hydrophones  had  a  much  smaller  array  baseline  than 
did  the  geophones,  calculations  implied  the  hydrophones  might 
serve  as  a  better  sensor  for  detection  and  measurement  of 
shallow  water  events,  at  least  for  ranges  closer  than  20 
kilometers.  This  conclusion  would  remain  valid  here,  even 
allowing  for  potential  improvements  in  geophone  SNR  expected 
to  be  achievable  with  proper  poliarization  filtering.  However, 
this  judgment  must  be  weighed  in  view  of  the  markedly  differ¬ 
ent  transfer  functions  of  the  two  types  of  sensors  and  of  the 
fact  that  both  sensor  signals  were  processed  at  the  same 
bandwidth.  It  should  also  be  cautioned  that,  owing  to  the 
relatively  small  hydrophone  array  baseline,  the  hydrophone- 
derived  SNR  and  correlation/beamforming  gain  are  subject  to 
large  uncertainties. 
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2.0  ANALYSIS 


The  investigation  discussed  in  this  section  was  conducted 
in  four  phases.  During  the  first  phase,  digital  representa¬ 
tions  of  the  shot  signals  were  graphically  analyzed  in  order 
to  develop  a  preliminary  description  of  the  seismic  propaga¬ 
tion  environment.  Four  different  shots  were  selected  for 
further  processing.  In  the  second  phase,  these  selected  shots 
were  subjected  to  time  series  and  spectral  analysis,  to 
determine  their  temporal  and  spectral  content.  In  the  third 
phase,  the  selected  signals  were  analyzed  by  means  of  coherent 
crosscorrelation  processing  to  enable  event  mode  detection 
and  isolation,  and  measurement  of  mode  propagation  direction 
(coherent  beamforming)  and  speed.  Estimates  were  obtained  for 
input  mode  SNR's  and  correlation  processing  gain.  Equivalent 
processing  was  then  performed  on  polarization  filtered  time 
series,  in  the  fourth  phase. 

2.1  Seismic  Modeling 

Low  resolution  digital  time  series  displays  of  the  geo¬ 
phone  and  hydrophone  signals  were  used  to  assemble  a  reduced 
travel-time  versus  distance  plot  for  all  the  shots  listed  in 
Table  2.  Truncated  versions  of  these  plots  are  presented  in 
Figures  2  through  7.  All  even  numbered  Figures  contain  sig¬ 
nals  only  from  shots  exploded  on  the  ocean  bottom.  Odd 
numbered  Figures  contain  signals  only  from  shots  exploded  at 
mid-depth  (about  20  m) . 

Seismic  mode  trajectories  were  identified,  wherever 
possible,  and  measured  from  the  time-distance  plots.  The 
Figures  show  these  trajectories,  superimposed  on  the  time 
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TABLE  2 

SHOT  DESCRIPTIONS 


1 


Shot 

■mill 

Location 

Shot 

Location 

A1 

108.89 

M 

A16 

14.96 

B 

A2 

109.04 

M 

A17 

11.17 

M 

A3 

109.13 

M 

A18 

11.24 

B 

A4 

101.96 

M 

A19 

9.28 

M 

AS 

92.38 

M 

A20 

9.27 

B 

A6 

83.36 

M 

A21 

8.85 

B 

A7 

73.67 

M 

A22 

8.26 

B 

A8 

64.89 

M 

A23 

7.79 

B 

A9 

55.85 

M 

A24 

7.37 

B 

A10 

46.37 

M 

A25 

6.93 

B 

All 

36.89 

M 

A26 

6.48 

B 

A12 

27.84 

M 

A27 

6.03 

B 

A13 

18.44 

M 

A28 

5.64 

B 

A14 

18.27 

B 

A30 

4.53 

B 

A15 

14.88 

M 

A32 

3.73 

B 

M  =  mid-depth  (^18m) 
B  =  bottom  (MOm) 

(From  Reference  1) 
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Figure  2  Reduced  Travel  Time - Di stance  Plot 
Geophone  3X  (Bottom  Shots) 
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Figure  3  Reduced  Travel  Time-Distance  Plot 
Geophone  3X  (Mid-depth  Shots) 
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series  displays.  A  two-layer  model  was  then  developed  to 
accommodate  the  observations.  Equations  necessary  to  compu¬ 
tation  of  the  seismic  model  parameters  are  discussed  below. 


Propagation  time,  t  ,  for  a  mode  which  travels  in  a 

ci 

uniform  layer  is 


t 


a 


(1) 


where  R  is  the  range  from  source  to  receiver.  Propagation 
time,  t  ^  ,  for  a  mode  which  traverses  from  the  top  of  the 
first  layer,  critically  refracts  into  a  second  layer,  then 
refracts  back  through  the  first  layer,  is 


where  h  is  the  first  layer  thickness;  c.  is  the  mode  speed 
in  the  first  layer;  c^  is  the  mode  speed  in.  the  second  layer. 
Arrival  time  difference,  At(a,afa),  between  the  modes  des¬ 
cribed  by  equations  (1)  and  (2)  ,  is 


At  (a,  afa)  =  R 

ca  cf 


(3) 


Arrival  time  difference,  At (afa,  aga),  for  two  modes  traversing 
as  described  by  equation  (2) ,  but  having  different  second 
layer  speeds,  c £  and  c  ,  is 

o 
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Finally,  the  arrival  time  difference,  At(afa,bfb),  for  the 
two  modes  traveling  as  described  by  equation  (2) ,  but  having 
different  first  layer  speeds,  ca  and  cb ,  is 


It  is  found,  from  equations  (3),  (4)  and  (5),  that 


6 (At (a ,afa) ) 
5R 


(3)' 


5 (At (afa ,aga)) 
6R 


(4)’ 


6(At  (afa  ,bfb) ) 
6R 


(5)’ 
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Equations  (3)  and  (4)  simplify  when  t  =  0,  a  condition  which 
can  be  established  on  time-distance  plots  by  projecting 
appropriate  mode  trajectories  to  their  intersection  point 
along  the  time  axis. 

Equations  (1)  through  (5)  and  their  derivatives  were 
used  to  compute  seismic  propagation  speeds  and  layer  thick¬ 
ness,  given  measurements  made  from  the  reduced  travel-time 
versus  distance  plots.  The  results  of  these  computations 
are  given  in  Table  3.  It  was  assumed  throughout  that  the 
hydroacoustic  propagation  speed  was  1.51  km/sec.,  as  reported 
in  reference  1,  from  earlier  experiments.  The  upper  layer 
shear  speed  was  deduced  indirectly  from  the  arrival  time 
difference  between  P-^2^1  and  •  Evidence  is  seen  for 

the  direct  arrival  of  ,  especially  in  the  z-axis  time 
traces,  but  it  is  not  as  consistent  or  obvious  as  for  the 
other  modes.  A  curious  effect  is  seen  wherein  the  first 
arrival  trajectory  for  shots  A32  and  A30  appears  delayed 
about  .1  seconds,  for  the  more  distant  shots.  This  is  as¬ 
cribed  to  a  low  velocity  zone,  situated  between  the  upper 
and  lower  layers,  which  must  have  been  encountered  for  ranges 
greater  than  about  5  kilometers.  Where  possible  in  the  cal¬ 
culations,  the  effects  of  the  low  velocity  zone  were  removed. 
Other,  much  lower  speed,  modes  were  perceived  in  the  nearer 
shots,  but  were  not  able  to  be  analyzed. 

High  resolution  time  series  displays,  shown  in  Figures  8 
through  23,  were  generated  for  signals  from  shots  A27,  A25, 
A20,  and  A14.  Note  that  the  waterwave  has  been  clipped  to 
maintain  a  substantial  amplitude  dynamic  range  for  the  seismic 
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TABLE  3 

DERIVED  SEISMIC  PROPAGATION  PARAMETERS 


l 

Propagation  Speed  fkm/sec 

Thickness  (km) 

Shear 

First 

.63  ±  .09 

1.74  ±  .02 

1 

1.32  ±  .12 

Second 

3.28  ±  .16 

2.34  ±  .06 

1.024  sec. 

Figure  8a  Shot  A27,  Geophone  5X 


1.024  sec. 

Figure  8b  Shot  A27,  Geophone  5Z 
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Figure  8c  Shot  A27,  Hydrophone  5 
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Figure  8e  Shot  A27,  Geophone  4Z 
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Figure  8f  Shot  A27,  Hydrophone 
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Figure  9a  Shot  A25,  Geophone  5X 
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Figure  9d  Shot  A25,  Geophone  4X 
-  27  - 
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Figure  9f  Shot  A25,  Hydrophone  4 
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Figure  10b  Shot  A20,  Geophone  5Z 
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Figure  lOd  Shot  A20,  Geophone  4X 
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Figure  11c  Shot  A14,  Hydrophone  5 
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Figure  lid  Shot  A14,  Geophone  4X 
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Figure  lie  Shot  A14,  Geophone  4Z 
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Figure  Ilf  Shot  A14,  Hydrophone  4 
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components.  In  each  of  the  signal  traces,  the  long  term 
mean  has  been  removed.  Transparencies  of  these  high  resolu¬ 
tion  traces  were  overlaid  to  facilitate  preliminary  judgments 
regarding  mode  polarization  states  and  array  transit  times. 

This  examination  of  polarization  led  to  the  assignment  of 
^^2^1  as  the  propagation  mode  for  the  third  seismic  arrival. 

It  also  was  helpful  in  searching  for  the  direct  S1  (?)  propa¬ 
gation  mode.  Seismic  energy  was  observed,  in  the  geophone 
signals,  to  persist  for  a  significant  time  after  the  water- 
wave  arrival.  Most  of  this  "coda"  appeared  to  arrive  as  a 
P  wave,  except  for  the  possible  S.  mode.  The  mechanisms  and 
propagation  paths  for  the  coda  were  not  pursued,  owing  to 
its  complexity.  A  listing  of  measured  signal-to-noise  ratios 
for  the  seismic  modes  is  presented  in  Table  4.  SNR's  were 
derived  using  maximum  signal  and  noise  amplitudes,  respectively. 

During  the  preliminary  investigation  of  the  high  resolu¬ 
tion  time  traces,  it  became  obvious  that  the  bulk  of  the  pre- 
first  arrival  signals,  in  any  given  shot,  were  very  highly 
correlated,  at  zero  relative  time  shift,  for  all  sensors  in 
the  array.  Correlated  amplitude  "glitches"  were  also  seen 
occurring  throughout  the  data  in  a  periodic  fashion  (about 
every  3/4  seconds) .  Since  a  signal  of  this  sort  could  arise 
naturally  only  as  a  plane  wave  incident  vertically  on  the 
array,  it  is  quite  unlikely  that  its  origin  lay  in  pheno¬ 
mena  typically  expected  to  characterize  the  array  environment. 
This  signal  is  discussed  further  in  later  sections. 

2 . 2  Spectral  Analysis 

The  selected  digital  data  signals  were  next  analyzed  to 
determine  their  temporal  and  spectral  content.  The  principal 
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TABLE  4 


MEASURED  SEISMIC  MODE  SIGNAL-TO-NOISE  LEVELS 

FOR 

VARIOUS  SENSORS 


_ _ Sensor* 

Shot  I 

Number  Mode  5X_ 5 Z_ 5H_ 4X_ 4 Z_ 4H 

A27  piP2Pl  4  3  5  1  ^  '6  6 

P1S2P1  17  ( 4 )  12  ( 2 )  16  (2)  10  (2)  3  (2)  16  (3) 

S1S2S1  15  (5)  9  (2)  14  (3)  12  (2)  1  (2)  14  (3) 


Pj  18  (4)  11  (2)  12  (3)  12  (2)  -3  (2)  12  (3) 


A25 

P1P2P1 

0  CD 

1  ( 2) 

7  (2) 

-4  (3) 

-13  (2) 

3  (3) 

P1S2P1 

17  (D 

5  (2) 

13  (2) 

3  (2) 

'1  (2) 

9  (2) 

S1S2S1 

17  (D 

4  (2) 

15  (2) 

2  (2) 

-3  (2) 

5  (2) 

P1 

27  (3) 

9  CD 

15  (3) 

13  (2) 

4  (2) 

13  (3) 

A20 

P1P2P1 

2  CD 

-5  CD 

-6  (D 

-13  CD 

-16  C  2 ) 

-12  (2) 

P1S2P1 

8  CD 

1  CD 

-2  (2) 

-4  C2) 

-14  (2) 

-l  CD 

S1S2S1 

8  (O 

6  C2) 

o  CD 

o  C  2) 

1 

O 

✓ — \ 
CO 

-l  CD 

P1 

21  (O 

11  CD 

3  CD 

5  C2) 

-6  (D 

2  CD 

A14 

P1P2P1 

2  CD 

-3  CD 

-14  CD 

-8  (2) 

-13  (2) 

P1S2P1 

3  CD 

-4  CD 

-10  (2) 

-6  CD 

1 

o 

/ - N 

-13  (2) 

S1S2S1 

5  CD 

-4  CD 

-10  CD 

-8  CD 

-7  CD 

-13  CD 

P1 

9  CD 

2  (2) 

3  CD 

-2  CD 

-4  (2) 

0  (2) 

*A11  signal-to-noise  levels  in  decibels;  estimated  error  in 
parenthesis . 
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working  tool  to  facilitate  this  process  was  the  spectrogram, 
which  is  a  graphic  representation  of  signal  power  as  a  func¬ 
tion  of  time  and  frequency.  An  extensive  collection  of 
spectrograms  for  the  selected  shot  signals  are  shown  in 
Figures  12  through  15.  Prior  to  spectral  analysis,  the 
signals  were  digitally  filtered  and  decimated  to  a  sampling 
rate  of  500  Hz. 

Every  display  line  in  each  spectrogram  was  developed 
from  1/4  second  segments  of  the  input  time  series.  Each  such 
data  block  was  augmented  with  zeros  to  a  total  length  of  one 
second.  Thus,  the  spectral  sampling  increment  is  one  hertz. 

The  temporal  analysis  window  (of  1/4  second  length)  was 
incremented  (slipped)  by  1/8  second,  thereby  producing  a  50% 
overlapped  spectrogram.  Analysis  bandwidth  spans  the  range 
from  -250  Hz  to  +250  Hz  (i.e.,  to  the  Nyquist  sampling  limits). 
After  the  Fourier  transform  of  each  zero-filled  time  data 
block  was  obtained,  its  logarithmic  power  spectrum  was  com¬ 
puted.  The  upper  10  dB  of  each  such  spectrum  was  converted 
into  a  10  level  density  format  for  plotting,  all  values  below 
(relative)  -10  dB  having  been  blanked.  Note  that  each  spec¬ 
trogram  line  is  therefore  scaled  only  locally,  not  globally. 
Also,  since  the  digital  signals  were  not  bandshifted,  the 
spectrograms  are  symmetric  about  0  hertz.  Indicators  on  the 
spectrograms  mark  the  onset  of  different  modes  which  comprise 
the  seismic  and  hydroacoustic  signal  structure. 

The  spectrograms  demonstrate  clearly  that  the  waterwave 
has  a  very  broadband  structure,  compared  with  the  seismic 
modes,  on  both  the  geophones  and  the  hydrophones.  Seismic 
energy  appears  to  be  confined  generally  to  less  than  25  Hz, 
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Figure  1 2u  Spectrogram,  Shot  A27,  llydroph 


I'igurc  12c  Spectrogram,  Shot  A27,  Gcophone  4Z 


Figure  1 2d  Spectrogram,  Shot  A27,  Hydrophone 
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Figure  1 2e  Spectrogram,  Shot  A27,  neophone  6X 


Figure  13a  Spectrogram,  Shot  A25,  Hydrophone 
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Figure  13b  Spectrogram,  Shot  A25,  Geophone  4X 


Figure  13c  Spectrogram,  Shot  A25,  Geophone  4Z 
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Figure  13d  Spectrogram,  Shot  A25,  Hydrophone 


Figure  13e  Spectrogram,  Shot  A25,  Geophone  6X 


Figure  14a  Spectrogram,  Shot  A20,  Geophone  5X 


Figure  14b  Spectrogram,  Shot  A20,  Geophone  SZ 


Spectrogram,  Shot  A20,  Hydrophone 


Figure  14d  Spectrogram,  Shot  A20,  Geophone  4X 


Figure  14e  Spectrogram,  Shot  A20,  Geophone  4Z 
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Figure  14f  Spectrogram,  Shot  A20,  Hydrophone 


Figure  14g  Spectrogram,  Shot  A20,  Geophone  6X 


Spectrogram,  Shot  A14,  Geophone  5X 
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Figure  15b  Spectrogram,  Shot  A14,  Geophone  5Z 
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Figure  15c  Spectrogram,  Shot  A14,  Hydropho 


Figure  15d  Spectrogram,  Shot  A14,  Geophone  4X 
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Figure  15e  Spectrogram,  Shot  A14,  Geophone  4Z 
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Figure  ISf  Spectrogram,  Shot  A14,  Hydrophone 
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Figure  15g  Spectrogram,  Shot  A14,  Geophone  6X 
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while  the  waterwave  energy  fills  the  analysis  bandwidth 
(limited  only  by  the  anti-alias  digital  filter  passband)  . 

A  moderate  degree  of  frequency  diversity  characterizes  the 
different  seismic  modes.  Signals  from  two  seismic  modes 
are  seen  to  occur  nearly  simultaneously,  and  to  overlap  the 
waterwave.  In  all  of  the  shots,  seismic  energy  can  be 
observed  for  a  significant  time  following  the  waterwave 
arrival.  Note  the  strong  5  Hz  noise  in  the  geophone  spectro¬ 
grams.  A  20  Hz  modulation  is  apparent  in  the  waterwave  power 
spectra,  implying  the  existence  of  a  time  delayed  arrival. 

The  approximate  time  delay  of  .05  seconds  (=  1/20  Hz)  corres¬ 
ponds  with  a  75  meter  difference  of  travel  path  length.  This 
could  be  attributed  to  a  shot  self - illumination  process,  in 
which  a  portion  of  the  hydroacoustic  signal  reflects  off  the 
ocean-air  interface  and  scatters  off  the  shot  bubble.  The 
multipath  could  also  have  been  caused  by  an  echo  near  the 
receiving  array.  However,  the  phase  and  spatial  consistency 
of  this  spectral  modulation  favor  the  former  explanation. 

Of  great  concern  to  this  signal  processing  investigation 
was  the  appearance  of  very  strong  harmonic  lines  in  all  of 
the  spectrograms.  At  first  glance,  the  line  structure  can 
be  seen  to  be  quite  similar  to  that  generated  typically  by 
marine  diesel  engines,  especially  those  having  a  small  number 
of  cylinders.  Upon  further  study,  it  was  found  that  this 
signal  occurred  simultaneously  and  highly  correlated  at  all 
sensors  throughout  the  data.  The  implication  was  that  such 
a  signal  had  moved  vertically,  in  a  plane  wavefront,  across 
the  array,  an  unlikely  circumstance.  A  significant  potential 
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exists  that  this  signal  was  extraneous  to  the  experimentally- 
acquired  signals,  having  its  cause  elsewhere.  The  strength, 
structure,  and  persistence  of  this  signal  also  had  adverse 
effects  on  subsequent  measurements. 

2 . 3  Coherent  Crosscorrelation  and  Beamforming 

Signals  received  at  the  array,  from  the  selected  shots, 
were  analyzed  using  coherent  crosscorrelation  techniques,  in 
order  to  detect  and  isolate  arrival  modes,  to  measure  their 
arrival  direction  and  apparent  array  transit  speed  (beam- 
forming)  ,  and  to  establish  input  SNR  and  output  processing 
gain.  A  discussion  of  the  mechanisms  by  which  the  coherent 
correlation  obtains  these  results  is  presented  first. 

The  coherent  crosscorrelation  function,  x>  of  two  digi¬ 
tal  time  series,  f  and  g,  is  defined  here  to  be 


X ( £ ,M,N)  = 


,  M+N-l  * 

N  £  fi  gi  +  J 
1  =m 


where 


M+N-l  *  4  M-H+N-l  * 

fj  fji  '£  *k  *k> 

j=m  J  J  k=M-£ 


Subscripts  in  equations  (6)  and  (7)  serve  to  label  elements 
of  the  two  digital  time  series.  Integration  time  of  the 
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correlation  is  equal  to  NSt,  where  6t  is  the  sampling  time 
increment.  Start  time  of  the  correlation  is  then  M6t,  with 
respect  to  the  first  series,  f.  The  relative  shift  (search) 
parameter  is  &,  corresponding  with  a  time  shift  of  £6t. 

If  the  two  signals,  f  and  g,  have  a  common  signal  com¬ 
ponent,  s,  then  the  time  series  can  be  redefined  by 

f  =  s+u;g  =  s+v,  (8) 

where  the  relative  time  shift  of  s  has  been  set  to  zero. 
Introducing  equation  (8)  into  equation  (6)  ,  one  has 


M+N-l 


X(St,M,N)  -  £  £  Cs. 

"  i=M 


i  +  2,' 


v .  „  +u  . 

i  +  £  l 


If  u  and  v  are  essentially  uncorrelated  with  s  at  the  time 
shift  which  matches  the  s  components  (i.e.,  at  zero  time 
delay  here)  ,  then  the  zero  lag  correlation  amplitude  is 
found  to  be,  from  equation  (9), 


1  M+N-l 

X ( 0 ,M ,N)  =  £  (s’?  s-  +  u*.  v.)  (10) 

"  i=M  11  11 

Further,  if  u  is  essentially  uncorrelated  with  v  at  zero  lag 
then  the  correlation  amplitude,  equation  (10),  has  the  value 

1  M+N-l 

X ( 0 ,M,N)  »  f  £  (s*  s.)  (11) 

"  i-M  1  1 
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where 


M+N-l  \  M+N-l  i 

*•[  Z  (s*Si+u*u.)]  [  £  (s£  sk+v*  v.)]  (12) 

j  =M  j  J  J  J  k=M  K  K  K  K 

Equation  (11)  can  be  rewritten  simply  as 

X (0 ,M,N)  =  {(1+SNRf-1)1  (1+SNRg-1) ^ }  (13) 

where  SNR^,  SNR  refer  to  the  signal-to-noise  ratio  of  the 
& 

first  and  second  time  series,  respectively.  If  both  channels 
have  identical  SNR,  then  equation  (13)  reduces  to 

X (0 ,M,N)  -  (1  +  SNR'1)'1 

or, 

SNR  =  (14) 

X 

by  which  the  SNR  can  be  measured.  Gain  of  the  correlation 
processing  is  defined  here  by  the  deflection  ratio,  D, 
where 

Z,  (15) 

ax 

In  equation  (15)  ,  x  is  the  mean  value  of  the  correlation 
function  (i.e.,  of  the  set  of  correlation  magnitudes  for 
many  relative  time  shifts)  and  is  the  standard  deviation 
of  the  correlation  function.  Equation  (15)  implies  that 
a  gain  is  achievable  only  if  the  common  signal  correlation 
magnitude  exceeds  the  mean  correlation  magnitude.  It  is 
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possible  to  compute,  from  a  theoretical  foundation,  expected 
values  for  the  mean  correlation  amplitude,  x>  and  the  correla¬ 
tion  standard  deviation,  a  ,  associated  with  a  pair  of 

A 

Gaussian  noise  input  signals.  These  values  can,  in  turn,  be 
used  to  quantify  the  maximum  achievable  correlation  process¬ 
ing  gain.  Previously,  it  was  shown  by  Carter  (Reference  2) 
that  the  magnitude  squared  correlation  probability  density 
function,  pN(u) ,  for  such  input  noise,  is 

PN(u)  =  (N-l)  (l-u)N"2  (16) 

where  N  is  the  number  of  (complex)  time  series  points  corre- 

lated,  and  u  *  x  •  For  convenience  in  notation,  it  has  been 

assumed  that  x  represents  the  modulus  (i.e.,  the  magnitude) 
correlation.  From  equation  (16)  ,  it  follows  that 

PN(X)  =  2(N-1)x(1-X2)N'2  (17) 

Consequently,  from  equation  (17), 

1 

XPN(x)dX 

*  \  /^r(N)/r(N+d  (18) 

If  Stirling's  expansion  for  the  gamma  function  is  utilized 
in  equation  (18)  ,  one  finds  that 
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Y  *V» 


(19) 


However,  the  expression  in  braces  approaches  1  //e,  for  large 
N,  so  that  equation  (19)  can  be  simplified  to  the  form 


x  %  W  5 

It  is  found  easily  that 

T  f1  2  ,  ,  .  1 

X  J  X  p(x)dx  =  JJ 


(20) 


(21) 


The  standard  deviation, 
(20)  and  (21) 


can  be  obtained  from  equations 


(22) 


Equations  (20)  and  (22)  can  be  used  to  determine  the  maximum 

deflection  ratio,  o  , 

’  m 


(23) 
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In  the  present  investigation,  correlation  lengths  of  64-200 
points  were  used,  so  that  the  expected  ranges  of  correlation 
output  parameters  were 

.06  <  x  <  .11 


.03  <  a  <  .06  (24) 

X 

24dB  <  D  <  29dB 
m 


Coherent  crosscorrelation  also  served  as  a  coherent 
beamformer  and  velocity  filter.  To  see  this,  consider  the 
geometric  array  configuration  shown  in  Figure  lb.  If  the 
incident  signal,  indicated  by  the  arrow  in  the  Figure, 
traverses  the  array  as  a  plane  wave,  then  the  arrival  angles 
are  related  to  the  arrival  time  delays  by  the  expression 


dijcoseij  .  ctjj 


Vri 


(25) 


In  equation  (23),  c  is  the  apparent  transit  speed,  is 

the  measured  relative  arrival  time  delay  as  observed  between 
two  sensors,  i  and  j;  and  r^  is  the  range  from  the  (apparent) 
source  to  sensor  k. 


Since 


013  =  612  ' 


then 


cosQ^j  =  cos<{icos012  +  sin<t>sin6^2 


(26) 
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If  equation  (25)  ,  for  i  =  1  and  j  =2,  is  introduced  into 
equation  (26),  one  finds  that 


sine 


12 


sinsj) 


rTl3 

t3n 


12 

12 


cos<{>] 


(27) 


Thus,  from  equations  (25)  and  (27), 


tane12 


t^) 

t12 


cos<f>] 


(  28) 


Equations  (25),  (27),  and  (28)  demonstrate  that  a  group 
of  coherently  related  time  delays  can  be  used  to  determine 
the  apparent  propagation  direction  and  speed  of  a  mode  which 
traverses  an  array.  It  is  observed,  from  equation  (25),  that 

Tij  *  Tjk =  Tik  (z« 

for  any  coherent  signal  which  might  traverse  the  array. 
Equation  (29)  links  multiple  pairwise  correlation  output 
data  together,  coherently.  Note  that  the  derivations  here 
apply,  except  of  course  where  <t>  =  180°  (linear  array),  to 
any  collection  of  sensors  selected  from  the  array. 

Two  types  of  display  were  utilized  to  present  the  raw 
correlation  output.  One,  the  density  plot  (Figures  16-19), 
is  analogous  to  the  spectrogram  density  plot.  Here,  the 
horizontal  axis  is  a  measure  of  relative  time  delay,  rather 
than  of  frequency.  Since  the  correlation  function  was  nor¬ 
malized,  its  display  was  globally,  rather  than  locally, 
scaled.  Also,  the  density  levels  were  derived  from  the 
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correlation  modulus,  rather  than  from  the  magnitude  modulus 
squared.  A  variety  of  A-line  correlation  format  displays 
are  included  in  Figures  20-26,  to  illustrate  detailed  behav¬ 
ior  of  the  coherent  correlation. 

Correlation  peaks  were  selected  from  the  data  for 
signals  arriving  from  the  direction  of  the  shots.  Correla¬ 
tion  peak  occurrence  times  were  used  to  identify  the 
respective  seismic  modes.  After  this  selection  process, 
the  apparent  propagation  speed  of  the  modes  were  calculated. 

SNR  and  correlation  gain  parameters  were  also  developed  from 
the  selected  events.  These  results  are  shown  in  Tables  5  and  6. 
The  bulk  of  the  seismic  parameters  reported  in  Tables  5  and  6. 
were  derived  from  x-axis  signals,  although  the  z-axis 
signals  yielded  similar  estimates  (as  is  evidenced  by  the 
z-axis  correlation  density  plots).  Regrettably,  the  geophone 
number  6  digitized  time  series  were  found  to  be  saturated 
at  maximum  negative  amplitude.  This  seriously  hindered  the 
independent  analysis  of  vertical  geophone  array  signals  in 
terms  of  beamforming  and  of  related  calculations. 

Hydrophone  results,  shown  in  the  Tables,  were  derived 
principally  from  hydrophone  sensors  number  4  and  5.  Since 
our  method  of  computing  the  digital  crosscorrelation  function 
had  a  precision  limited  by  the  sample  time  increment,  it  was 
found  that  estimates  for  propagation  direction  and  speed, 
based  on  the  small  hydrophone  array,  were  relatively  coarse 
compared  with  the  geophone-based  estimates.  For  this  reason, 
the  hydrophone  measurements  were  comparatively  less  accurate. 

All  of  the  results  in  Tables  5  and  6  represent  the  best  data 
obtained  using  several  different  sets  of  correlation  process¬ 
ing  parameters. 
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TABLE  5 

CORRELATION  PROCESSING  GAIN  AND  SIGNAL-TO-NOISE  LEVELS 


Deflection  Ratio  (dB) 

Signal- 

to-Noise  Ratio  (dp) 

Shot 

Number 

Mode 

Ml 

5X/4X 

5X/6X 

5H/4H 

A27 

P1P2P1 

11, 

6 

9, 

i 

12, 

13 

8, 

17 

P1S2P1 

11, 

14 

6, 

is 

6, 

13 

9, 

38 

S1S2S1 

11, 

14 

5, 

16 

6, 

13 

9, 

38 

P1 

10, 

2 

- 

- 

1, 

0 

6, 

34 

Water 

18, 

10 

17, 

8 

18, 

12 

21, 

35 

A25 

P1P2P1 

10, 

0 

9, 

10 

11, 

13 

9, 

10 

P1S2P1 

10, 

9 

9, 

6 

9, 

20 

8, 

26 

S1S2S1 

6, 

9 

3, 

1 

9, 

15 

9, 

30 

P1 

8, 

14 

- 

' 

- 

- 

10, 

31 

Water 

18, 

12 

15, 

6 

.  .  17. 

14 

18, 

44 

A20 

P1P2P1 

- 

- 

- 

- 

- 

- 

- 

- 

P1S2P1 

12, 

5 

10, 

7 

8, 

10 

10, 

6 

S1S2S1 

U, 

6 

8, 

7 

8, 

8 

7, 

20 

P1 

12, 

10 

10, 

9 

10, 

6 

10, 

14 

Water 

18, 

8 

17, 

S 

17, 

11 

20, 

38 

A14 

P1P2P1 

- 

- 

- 

1, 

-1 

- 

- 

P1S2P1 

8, 

1 

- 

- 

- 

- 

- 

S1S2S1 

2, 

5 

- 

- 

- 

- 

- 

P1 

10, 

8 

- 

- 

- 

6, 

20 

18, 

3 

2 

20, 

14 

19, 

34 
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Correlation  Density  Display, 
Geopnones  5X  vs.  6X 
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Figure  16c  Correlation  Density  Display,  Shot  A27 , 
Gcopiiones  4Z  vs.  5Z 
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Figure  I6d  Correlation  Density  Display,  Shot  A27, 
Hydrophones  4  vs.  5 
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Figure  17a  Correlation  Density  Display,  Shot  A25, 
Geophones  4X  vs.  5X 
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Figure  17b  Correlation  Density  Display,  Shot  A25, 
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Figure  17e  Correlation  Density  Display,  Shot  A25, 
Hydrophone  4  vs.  Geophone  5X 
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Geophoncs  5X  vs.  6 
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Figure  18c  Correlation  Density  Display,  Shot  A20, 
Geophonos  M  vs.  5Z 


Figure  18d  Correlation  Density  Display,  Shot  A20, 
Hydrophones  4  vs.  5 
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Figure  18e  Correlation  Density  Display,  Shot  A20, 
Hydrophone  4  vs.  Geophone  5X 
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Figure  19c  Correlation  Density  Display,  Shot  A14 
Gcoplioncs  4Z  vs.  5Z 
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Figure  19d  Correlation  Density  Display,  Shot  A14, 
Hydrophone  4  vs .  5 
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Figure  21a  Correlation  Function,  Shot  A25,  Geophones  5X  vs.  6X 


Figure  21c  Correlation  Function,  Shot  A25,  Hydrophon 
Mode  P,S0P.. 


Correlation  Function,  Shot  A20,  Geophones  5X  vs.  4X 
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igure  23b  Correlation  Function,  Shot  A27,  Geophones  5X  vs.  6X 
Water  Wave 
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Figure  24b  Correlation  Function,  Shot  A25,  Geophones  5X  vs.  6X, 
Water  Wave 
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Figure  25b  Correlation  Function,  Shot  A20,  Geophones  5X  vs 
Water  Wave 


gure  25c  Correlation  Function,  Shot  A20  ,  Geophones  5Z  vs.  4Z 
Water  Wave 


tfl 

> 

X 

LO 

to 

<D 

C 

o 

X 

a 

o 

<u 

a 


h- 

rH 

< 


o 

to 


o 

G 

3 

tu 

c 

o  a> 

•H  > 
+J  CO 
cC  SC 

pH 

CD  G 
J-  (D 
G  4-> 
O  CO 

CJ  3: 


CO 

CM 

CD 

G 

3 

oo 


Figure  26b  Correlation  Function,  Shot  A14 ,  Geophones  SX  vs.  6X, 
Water  Wave 


Figure  26c  Correlation  Function,  Shot  A14,  Geophones  5Z  vs.  4Z, 
Water  Wave 
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As  expected  from  earlier  analysis,  the  correlations 
were  characterized  by  a  strong  peak  at  zero  relative  time 
delay,  for  all  sensor  pairs.  Related  strong  satellite 
correlations  appeared  at  intervals  of  about  .025  seconds. 
Weaker  correlations  can  be  seen  halfway  between  the  major 
peaks.  These  periodic  correlation  artifacts  are  readily 
attributable  to  the  interference  signal  multiple  tonal 
structure.  Assume  that  an  interference  signal,  f  (t),  consists 
of  at  least  two  tones,  having  frequencies,  and  cd2 ,  res¬ 
pectively.  That  is. 


iu),T  loj^T 

f(t)  =  e  1  +  e  1  (30) 

where  a  could  be  complex-valued.  The  long  term  crosscorrela¬ 
tion  function,  x("0»  obtained  from  this  signal  being  present, 
at  zero  relative  time  lag,  in  both  channels,  would  be 


x(t) 


(31) 


The  modulus  correlation  amplitude,  from  equation  (31),  is 


Ix(t)|  =  /  — - — y 
\ 1+ I a i 

which  indicates  that  the  correlation  function  should  be 
modulated  at  a  frequency,  Since  the  spectrograms 

show  the  presence  of  39  Hz  and  78  Hz  harmonic  component  in 
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d+ 
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1  +  2 
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'cos  ( 


Wo  -to 


)  t 


(32) 
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the  interference,  then  the  resultant  correlation  function 
should  have  been  modulated  at  a  period  of  1/39  Hz  =  .026 
seconds.  Presence  of  the  next  higher  harmonic,  also  seen 
in  the  spectrograms,  should  have  resulted  in  a  modulation 
period  of  1/78  Hz  =  .013  seconds.  Clearly,  the  observed 
strong  modulation  pattern  in  the  correlation  functions  were 
caused  by  the  interference  signal. 

One  explanation  for  such  an  interference  signal  is  that 
it  was  due  to  a  plane  wave  incident  simultaneously  across 
the  array;  that  is,  due  to  a  wavefront  arriving  from  either 
directly  above  or  below  the  array.  Since  this  signal  was 
present  for  all  shots  (covering  a  period  of  hours) ,  it  seems 
more  likely  that  it  was  caused  otherwise,  perhaps  during 
the  recording,  playback,  or  digitizing  stages.  Intriguingly , 
the  spectral  structure  of  the  interference  is  not  inconsis¬ 
tent  with  that  observed  from  diesel  marine  engines. 

In  essentially  all  results,  except  for  the  waterwave 
arrivals,  the  interference  signal  has  significantly  altered 
correlation  time  delays,  amplitudes,  and  statistics.  Conse¬ 
quently,  the  non-waterwave  data  in  Tables  5  and  6  should  be 
regarded  with  great  caution.  Consideration  was  initially 
given  to  filtering  out  the  interference  by  additional  digi¬ 
tal  processing  (e.g.,  bandlimiting ,  spectral  "notching," 
spectral  whitening) ,  but  it  was  recognized  that  this  would 
have  tended  to  modify  signal  properties  in  ways  which  could 
not  be  easily  gauged.  The  effort  required  to  quantify  these 
effects  was  too  extensive  to  undertake  within  the  limits  of 
this  investigation. 
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It  is  evident,  from  Table  5,  that  the  correlation- 
derived  SNR's  were  nearly  always  higher  for  the  hydrophones 
than  for  the  geophones.  This  probably  was  the  result  of  the 
fact  that  the  hydrophones  were  well  isolated  from  the  low 
frequency  (^5  Hz)  seismic  noise  which  affected  the  geophones. 
The  spectrograms  demonstrate  this  effect.  On  the  other  hand, 
correlation  gain  (deflection  ratio)  performance  was  found  to 
be  roughly  comparable  for  both  types  of  sensor.  This  was 
not  an  unexpected  result,  since  directional  properties  of 
the  arriving  signals  were  maintained  in  the  geophone  output, 
but  lost  in  the  hydrophone  output.  That  is,  in  contrast 
with  a  three-axis  geophone,  a  hydrophone  would  be  sensitive 
to  signals  (including  noise)  propagating  from  all  directions. 

No  clear  pattern  emerges  from  the  correlation-measured 
SNR  and  gain  parameters  by  which  one  could  make  decisive 
judgments  regarding  the  relative  merits  of  hydrophones  versus 
geophones  in  this  application.  Indeed,  the  directly-measured 
(Table  4)  and  correlation-measured  (Table  5)  SNR's  do  not 
reflect  usual  signal  propagation  attenuation  losses.  Note, 
however,  that  the  interference  signals  might  have  strongly 
masked  this,  and  perhaps,  other  properties  of  the  shot-caused 
experimental  signals. 

Apparent  array  transit  speeds  for  the  seismic  modes  are 
seen  to  be  larger  than  derived  from  the  time  series  analysis. 
For  the  purposes  of  the  seismic  modeling  analysis,  it  had 
been  assumed  that  all  signals  propagated  ideally  in  a  hori¬ 
zontal  layered  medium,  undergoing  transitions  at  critical 
angles.  Commonly,  apparent  seismic  transit  speeds  over  an 
array  do  not  match  actual  propagation  speeds,  owing  to  local 
layer  undulations  and  inhomogeneities.  Also,  a  non-critically 
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refracted  wavefront  incident  on  an  array  exhibits  a  larger 
apparent,  compared  with  real,  propagation  speed.  No  attempt 
was  made  to  model  potential  local  geologies  to  account  for 
the  apparent  speed  anomalies.  However,  the  waterwave  also 
yielded  an  apparent  speed  somewhat  larger  than  expected, 
although  it  was  not  subject  to  the  same  propagation  vagaries 
as  were  the  seismic  waves. 

2 . 4  Polarization  Filtering 

The  concept  of  polarization  filtering  is  based  on  the 
fact  that  a  seismic  signal,  as  observed  at  a  tri-axial  geo¬ 
phone,  must  generate  displacements,  along  the  axes,  that  have 
specific  relative  phase  relationships.  Consider  the  seismic 
signals  as  defined  in  Figure  2S .  If  a  compressive  wave  is 
incident  on  the  geophone  from  the  fourth  quadrant,  as  in 
diagram  (a),  the  x-axis  displacement  will  be  at  maximum  posi¬ 
tive  when  the  z-axis  displacement  is  at  maximum  negative. 

Thus,  a  superposition  of  the  x  and  z  time  traces  for  this 
signal  would  reveal  (in  the  absence  of  noise)  that  they  were 
out  of  phase  by  180°  (i.e.,  by  a  minus  sign).  A  shear  wave, 
incident  from  the  fourth  quadrant,  as  in  (b) ,  would  reveal 
x  and  z  time  traces  that  were  exactly  in  phase  (i.e.,  identi¬ 
cal  except  for  amplitude  scale).  A  compressive  wave,  incident 
from  the  first  quadrant,  as  in  (c),  would  yield  x  and  z  time 
traces  in  phase,  equivalent  to  that  in  the  fourth  quadrant 
shear  case,  (b) .  Rayleigh  modes  would  have  a  90°  or  270° 
phase  difference,  depending  on  their  propagation  direction 
and  type  (prograde  or  retrograde) .  Similar  arguments  can  be 
advanced  for  Love  and  other  waves.  Consequently,  processing 
gains  should  be  realizable  in  extracting  such  signals  from 
noise,  the  latter  by  definition  having  random  relative  phase 
among  the  three  sensor  measurement  axes. 
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The  polarization  filtering  techniques  used  here  operated 
in  the  signal  spectral  domain.  Fast  Fourier  transformation 
was  first  performed  on  the  three  axis  signals,  followed  by 
computation  of  their  relative  phase  differences  at  each  spec¬ 
tral  element  (picket).  Weights  were  applied  to  the  separate 
spectra  in  a  manner  which  favored  the  desired  phase  angles. 
Spectral  pickets  which  deviated  from  the  desired  relative 
phase  relationship  were  reduced  in  amplitude,  as  a  function 
of  their  deviation.  Polarization-filtered  time  series  were 
then  constructed  via  recombination  of  these  (time-overlapped) 
modified  spectra.  Fourier  transformation  integration  lengths 
were  chosen  in  such  a  way  as  to  provide  near-optimal  frequency 
and  phase  resolution,  for  the  various  arriving  modes.  Approxi¬ 
mately,  this  resolution  varies  as 


(N6t  /  /SNR) "1 
for  frequency  and 


2tt  /  (/SNR)  x 

for  phase,  where  N6t  is  the  integration  time  of  and  SNR  is 
the  signal-to-noise  ratio  in  the  analysis  time  window. 

Polarized  time  series  were  then  analyzed  in  the  same 
manner  as  had  been  the  unfiltered  time  series.  Spectrograms 
of  selected  filtered  time  series  are  presented  in  the  Figures, 
for  comparison  with  untiltered  time  series.  A  collection  of 
crosscorrelation  density  plots,  computed  from  polarization 
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Figure  28a  0°  Polarized  vs.  Unpolarized  Spectrogram, 

Shot  A27 ,  Geophone  5X 
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Figure  28b  180°  Polarized  vs.  Unpolarized  Spectrogram, 
Shot  A27 ,  Geophone  5X 


Figure  28c  0°  Polarized  vs.  Unpolarized  Spectrogram, 

Shot  A27 ,  Geophone  5Z 


Figure  28d  180®  Polarized  vs.  Unpolarized  Spectrogram, 

Shot  A27,  Geophone  5Z 


Figure  29a  0°  Polarized  vs.  Unpolarized  Spectrogram, 

Shot  A25,  Geophone  5X 
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Figure  29b  180°  Polarized  vs.  Unpolarized  Spectrogram, 

Shot  A25,  Geophone  5X 
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Figure  29c  0°  Polarized  vs.  Unpolarized  Spectrogram, 

Shot  A25 ,  Geophone  5Z 


Figure  29d  180°  Polarized  vs.  Unpolarized  Spectrogram, 

Shot  A25,  Geophone  5Z 
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Shot  A25,  Geophones  4Z  vs. 


filtered  data,  is  also  presented  in  the  Figures.  Unfortu¬ 
nately,  in  the  large  part,  the  polarization  processing 
performed  here  did  not  perform  well  at  all.  These  results 
were  in  spite  of  the  fact  that  well-defined  phase  relation¬ 
ships  can  be  seen  clearly  in  the  high  resolution  time  traces 
Figures  8  through  11.  It  is  possible  that  the  transfer  and 
modification  of  existing  polarization  filter  computer  code 
from  previous  applications  were  not  achieved  without  error. 

A  modest  degree  of  successful  polarization  filtering, 
nonetheless,  is  evident  in  spectral  and  correlation  displays 
In  the  cases  where  inphase  time  series  were  favored,  for 
the  x  and  z  axes,  the  hydroacoustic  signals  (arriving  as  a 
compressive  wave  from  the  first  quadrant)  were  preserved 
(note  the  persistence  of  the  inphase  interference  signal). 

On  the  other  hand,  where  180°  phase  separations  were  favored 
the  hydroacoustic  (and  interference)  signals  were  removed. 

No  obvious  pattern  emerged,  however,  for  the  compressive  and 
shear  seismic  modes,  by  which  to  permit  their  identification 
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